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A computational study was conducted to evaluate the effectiveness of using a viscous flow 
solution in an ice accretion code and the resulting accuracy of aerodynamic performance 
prediction. Ice shapes were obtained for one single-element and one multi-element airfoil 
using both potential flow and Navier-Stokes flowlields in the LEWICE ice accretion code. 
Aerodynamics were then calculated using a Navier-Stokes flow solver. 

Nomenclature 

two-dimensional 
section drag coefficient 
computational fluid dynamics 
section lift coefficient 
pressure coefficient 
liquid water content 
median volume diameter 
Mach number 

Reynolds-Averaged Navier Stokes 
Reynolds number 
Spalart-Almaras 
freestream velocity 
angle of attack 


I. Introduction 

A CCRETION prediction codes are improving their accuracy and productivity, yet some recognized 
shortcomings remain. Flowfields where viscous and unsteady effects are important may not be modeled 
correctly, such as cases where large regions of separated flow are present. These situations can occur in a variety of 
scenarios likely to contain the critical ice shape- large glaze ice horns, high angles of attack, and also in multi- 
element airfoils. With airfoil performance as a criteria, the trade between accuracy and effort still remains a 
question. Attempts to determine the usefulness of Navier-Stokes flow solutions in ice accretion prediction have only 
recently begun. This paper begins to address whether the additional effort of implementing a computational fluid 
dynamics (CFD) flow solver solution into the ice accretion process is truly significant in terms of performance 
degradation. 

While NASA has integrated a Navier-Stokes flow solver with its primary ice accretion code, LEWICE 1 , 
validation to date is somewhat limited, and results appear mixed. The question of whether the improvements in the 
predicted ice shapes warrant the additional computational expense has not been fully addressed. Further study is 
needed to determine if these improvements are significant in terms of performance degradation. 

With plans to incorporate a Navier-Stokes flow solver more directly into the ice accretion and droplet trajectory 
modules of an ice accretion code, it is important to quantify the gain in performance degradation accuracy that can 
be obtained. This will help NASA decide how tightly or loosely coupled the future generation of codes should be 
with the thermodynamic, trajectory and multi-phase elements of ice accretion prediction. 
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II. Single Element Airfoil 

A classical symmetric airfoil was chosen for the initial single-element phase of this study. A NACA 0012 airfoil 
section at 3.5 degrees angle of attack and a flight speed of 102.8 m/sec was used. The chord length was 0.5334m. 
These conditions correspond to a freestream Mach number of 0.324 and a Reynold’s number of 2,400,000. The 
static temperature was 250. 37K and the static pressure was 100,000 N/nr. The icing conditions were selected from 
the LEWICE validation database, and correspond to a seven-minute exposure in rime ice conditions. A median 
volume diameter of 2 Op and a liquid water content of 0.55 g/m 3 were chosen. 

These conditions correspond to run 405 of the 1998 Ice Shape Comparison test. This test to determine 
experimental ice shapes was conducted in the NASA Glenn Icing Research Tunnel using a model NACA0012 with 
a 21 -inch chord. The ice tracing used for comparison was obtained at a spanwise station of 36-inches. 

ICEG2D, the NASA Glenn integrated Navier-Stokes flow solver and ice-accretion code 2 , was then used in batch 
mode to calculate a 7-minute ice shape. ICEG2D initially distributes points along a user-defined airfoil geometry 
using a specified number of control points. These points are then automatically clustered in regions of high 
curvature to ensure geometric fidelity. A NURBS approximation of the iced airfoil is then generated using 
parameters suitable for either a smooth or rough surface, as defined by the user. A marching algorithm with a 
parabolic solver is used to generate a structured grid for each icing time step (typically 1 minute). 

A single-zone C-type grid with wake cut was developed. A downstream distance of 14 chord lengths was 
selected for the downstream zone boundary. A non-dimensional distance to the first grid point off the airfoil surface 
of 0.00001 was used to give a y+ value of approximately 1. Grid quality for the default input values in ICEG2D has 
already been established’ for typical cases. 

ICEG2D was run for this case in Level 0 coupling mode. In Level 0 coupling mode, the Navier-Stokes solution 
is used in place of the panel method solution in LEWICE. This affects both the droplet trajectories and the boundary 
layer computations in the ice accretion process. ICEG2D automatically supplies the flow solution and grid file 
needed by LEWICE and sets the appropriate input parameter (IGRID=1) when Level 0 is specified. The WIND code 
also provides surface heat transfer information, which can also be included in the ice accretion step through the 
LEWICE input parameter IHTC. This is known as Level 1 coupling, and also includes the Level 0 process. For the 
rime shape used in this study, using the surface heat transfer data from WIND did not appear to provide significant 
improvements over Level 0. This is not necessarily the case for all ice shapes. 

The 7-minute rime ice shape was then compared with the corresponding ice shape developed experimentally in 
the IRT and with a LEWICE-derived ice shape using potential flow. The ice shapes calculated by ICEG2D are 
shown at one minute intervals in Figure 1. 

This demonstrates accurate prediction of ice 
shape for cold conditions such as this. 

Aerodynamic performance was then 
calculated for the ice shapes derived from 
Navier-Stokes, potential flow and 
experiment. A flowfield for the clean airfoil 
was also computed. 

Aerodynamic performance was calculated 
using WIND 4 , a general-purpose flow solver. 

WIND is a second-order accurate finite- 
volume code which solves the Navier-Stokes 
equations in conservative form. The code has 
Runge-Kutta and Global Newton schemes for 
time accurate computations, but for this study 
the Reynolds Averaged Navier Stokes 
(RANS) scheme was used. 

In addition to the flow solver itself 
WIND includes a number of associated utilities, and several of these were used. The program cfcnvt converts 
between the standard PLOT3D format and the cgd format that is specific to WIND. The cfsubset allows a user to 
manipulate an existing grid by cropping various indices. The utility code gman is used to apply boundary conditions 
and zone connectivity information to a generic grid prior to input into WIND. 
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WIND allows the user to select from a 
variety of algebraic, one-equation and two- 
equation turbulence models. The Spalart- 
Almaras (S-A) one-equation turbulence 
model ’ was selected for its robust and general 
applicability to icing problems 6 . Some initial 
calculations were also made using the K-e 
model of Chien 7 , but no significant 
improvements were found for this case, so the 
S-A model was used throughout. 

The convergence history of the 7-minute 
ice shape is shown in Figure 2. The number of 
iterations recommended for use in ICEG2D is 
2000, but for this study solutions were run 
through 800 cycles of 20 iterations each. Run 
times to 16,000 iterations on a single- 
processor SGI were typically about 4 l A hours 
of CPU time. Convergence of the log(L2) 

residual was typically obtained to the fifth order in significantly less than 16,000 iterations. 

Table 1 shows the lift, drag and grid size for the clean, LEWICE, experimental and ICEG2D (taken at 1 -minute 
intervals) solutions. 



Figure 2. Convergence of NACA0012 Solution 
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347x101 
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Table 1. Section Coefficients Calculated for Single-element Airfoil Performance. 


Figures 3,4,5, and 6 show the total velocity field for the four cases. Of particular interest is the roughness in the 
experimental data and the corresponding effect this has on the section drag coefficient. The surface pressures were 
also calculated, and these are shown in Figure 7 for the three ice shape methods, along with the clean airfoil surface. 
Skin friction coefficient for the four cases is shown in Figure 8. 

It is now possible to run multiple time steps for single-element airfoils, and suitable results can be obtained for 
rime ice shapes at low to moderate angles of attack. But the subject still has not been completely investigated, and 
this paper represents only the preliminary phase of an investigation into the practical uses and applications of 
Navier-Stokes-coupled ice accretion prediction. The effects of wall temperature boundary condition (Level 1 
coupling) will be further investigated, for warm glaze ice shapes in particular. The flow solver itself may suffer from 
limitations at high angles of attack where available turbulence models don’t really handle stall very well 6 9 More 
advanced methods, such as DES and LES, need to be investigated further. The benefits of the added coupling may 
eventually play a role in the development of new or modified roughness models in the boundary layer computations. 

TIT. Multi-Element Airfoil 

Although a rime ice shape on a single-element airfoil can be a challenging aerodynamic problem to solve, it is 
not one which necessarily requires a viscous flow solution to obtain satisfactory results. Those ice shapes which are 
more dominated by viscous effects either exhibit large regions of flow separation caused by massive glaze-ice horns, 
or occur in airfoils which naturally have regions of reversed or separated flow, such as might occur in multi-element 
airfoils. 

A modern multi-element airfoil in the landing configuration was chosen for the second part of this study. The 
Douglas/NASA Lewis/NASA Langley Advanced High-Lift Multi-Element Airfoil Program 10 ' 11 was established in 
1993 to improve the understanding of ice accretion effects on multi-element airfoils. 

At the time these calculations were made, the multiple time-step ice accretion used in the single-element study 
was not available for multiple-element airfoils, so the ice accretion was limited to a single time step. However, 
LEWICE was capable of reading in multi-block flowfield solutions, so an overset chimera grid approach was taken. 
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Figure 3. Total Velocity (LEWICE-derived shape) Figure 4. Total Velocity (ICEG2D-derived shape) 




Figure 5. Total Velocity (Experimental Shape) 


Figure 6. Total Velocity (Clean Airfoil) 



Figure 7. NACA 0012 Pressure Coefficient 
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Figure 8. NACA0012 Skin Friction Coefficient 
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Conditions corresponding to Case 1 1 (a multi-element case) in the LEWICE user’s manual were selected, with 
one exception: a monodisperse droplet distribution was used as opposed to a Langmuir distribution. Otherwise, the 
conditions were the same. Freestream velocity was 88.5m/sec, corresponding to a Mach number of 0.27. MVD was 
2 Op and LWC was 0.60g/m 3 . Angle of attack was 8.0 degrees. The icing exposure time was 6 minutes at a 
freestream temperature of 268.15K. 

The process by which the multi-element ice shape was calculated consisted of the following steps: create x.y 
geometry files for the slat, main element and flap; generate 2D grids for slat, main element and flap using ICEG2D 
tools; convert grid ( *.grd ) files to WIND ( *.cgd) format using the WIND utility cfcnvt; modify slat and flap cgd files 
using the WIND utility cfsubset; combine slat, main and flap grids into one cgd file using cfcnvt; create 
holes/fringes/overlaps and assign boundary conditions using gman; scale the grid; run the flow solver; monitor 
convergence; post-process the solution; and enter the flow solution into the 2D version of LEWICE. 

The geometries for the slat, main element and flap were entered separately into the point-distribution and grid- 
generation modules using ICEG2D in interactive mode. Keeping the overall topology in mind, different values were 
used to create the slat and flap grids, as shown in Table 2. This produced essentially three different single-element 
airfoils which then had to be combined into a single grid. The 2D grid for the main element was used as the 
background grid. 

Conversion of the ICEG2D grid files into cgd format made use of the existing WIND utility, cfcnvt. The double- 
precision grid files were in Plot3D ASCII format, with the multi-zone, 2D and I-BLANK flags active. This step was 
necessary for all three elements. 



Slat 

Main element 

Flap 

distOuterB oundary 

0.5 

15.0 

0.5 

nLayers 

100 

100 

100 

distWake nl 

0.5 

14.0 

0.5 


Table 2. Grid Generation Inputs for Multi-element Airfoil. 


The utility cfsubset was then used to trim the outer I and J grid lines in the flap and slat grids. The gen grid 
command in ICEG2D was found to have difficulties converging if the grid boundaries were set too close to the 
airfoil. The flap grid should not overlap the slat grid at all, and vice-versa, so cropping was required. The flap and 
slat zones were trimmed so that both were totally enclosed by the main grid, yet neither intersected the main body by 
more than a few grid points in the interior (if any). It was also necessary not to trim so much that they didn’t have 
some overlap between the intended holes (ie. the fringe) and the outer zone boundary. Some engineering judgement 
was involved in selecting the I and J indices 
which were to remain after cropping. These 
indices were related to but not the same as 
the cutting surfaces that would used in gman 
in a later step. A J index of approximately 

0. 25 to 0.50 chord lengths and I indices 
approximately 5-10 grid points aft of the 
trailing edge were found to be sufficient. 

This step involved the most user intervention 
and would likely be the most difficult to 
automate, but was necessary to prevent 
undefined points in coupled boundaries in the 
next step. 

The slat, main and flap grids were then 
combined into a single grid in WIND cgd 
format using the cfcnvt utility. By 
convention, the main element was set as zone 

1, the slat as zone 2 and the flap as zone 3. 

The grid is shown in Figure 9. 

The most involved step was the creation 
of fringes, overlaps and boundary conditions 
in the gman utility. The process used for the 
three-element airfoil was similar to that used 
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Figure 9. Multi-Element Airfoil Grid 
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for an airfoil with flap in the NPARC Validation Archive 12 . The selection of cutting surfaces and subareas for 
boundary conditions was practically determined iteratively. The I indices to be used as cutting surfaces were within 
one grid point of the trailing edge or else on the trailing edge. The J index used as a cutting surface must be 
somewhat less than the maximum J index remaining after cropping. 

The automatic grid generator in ICEG2D worked best when the airfoil geometry was given in terms of the unit 
chord, in other words the trailing edge was at x=l in the x,y geometry file (*.c/«). It will work otherwise, for 
example, using wind-tunnel or full-scale geometries. However, the default values are based on a chord-normalized 
geometry, which will result in an improper Reynold’s Number in the flow solver unless the grid is scaled. In order to 
match Reynold’s number, the grid was scaled up from chord units to physical units (inches) using the global scale 
option in gmcin. Because LEWICE operates in chord-normalized dimensions, a later conversion would become 
necessary. 

Since only the point-distribution and grid-generation modules of ICEG2D were used in the multi-element case, a 
WIND input file had to be created. A grid of 50,000 points took approximately 8-10 hours to reach a converged 
solution on an Origin SGI. This corresponds to roughly 16-20K iterations. 

Post-processing the solution first required a conversion from c/7 to PLOT3D format, using the cfpost utility. A 
second program, a LEWICE utility specific to multi-element airfoils, was then required to convert the standard 
PLOT3D solution to LEWICE format. In addition, this utility performed several other functions. The flow solution 
had to be converted to appropriate physical units and the grid was also converted back to chord units. The three 
zones were re-ordered to reflect LEWICE input requirements. This convention required that the element geometry 
files be entered into LEWICE in the following order: slat, then main element, then flap. The viscous solution and 
grid files must also be renamed in order for LEWICE to accept them as inputs. Although it was not necessary in this 
case, binary to ASCII conversion is necessary if transferring the viscous WIND solutions to a PC rather than an SGI 
for running LEWICE. 

For multi-element geometries using viscous flowfield solutions, LEWICE cases currently only accept single 
time-step and single drop-size distributions. LEWICE still requires a basic geometry file to run, even when using a 
viscous solution. The geometric input files 
should correspond exactly to the surface 
definition in the grid files, ie. no offsets or 
smoothing to approximate surface streamlines 
should be made. 

The ice shapes predicted for the multi- 
element slat from both viscous and potential 
flow are shown in Figure 10, along with the 
experimental data. The predictions did not 
appear to be completely reliable, and thus 
further investigation is required. The 
discrepancy here is almost certainly due to the 
monodispersed droplet distribution that was 
used. Even at zero degrees angle of attack, the 
slat element (zone 3) did not converge as well 
as the upstream elements in zones 1 and 2. 

There may have been some unsteadiness 
involved, and improvements to the overall zone 
topology are likely required. 

It now appears feasible to automate the LEWICE suite of codes to provide for multiple time-step ice accretions 
on multi-element airfoils. However, additional improvements to the suite of computational tools are required before 
the process can be automated. These improvements include more robust tools for crossing between computing 
platforms, a multi-element capability for the automated grid generator, and full implementation of all LEWICE 
drop-size and droplet distribution capabilities in the viscous flow solution mode. 
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FigurelO. Multi-Element Ice Shape (Slat) 
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